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Abstract 



The shear viscosity for a heated granular binary mixture of smooth hard 
spheres at low-density is analyzed. The mixture is heated by the action of 
an external driving force (Gaussian thermostat) which exactly compensate for 
cooling effects associated with the dissipation of collisions. The study is made 
from the Boltzmann kinetic theory, which is solved by using two complemen- 
tary approaches. First, a normal solution of the Boltzmann equation via the 
Chapman-Enskog method is obtained up to first order in the spatial gradients. 
The mass, heat, and momentum fluxes are determined and the corresponding 
transport coefficients identified. As in the free cooling case [V. Garzo and 
J. W. Dufty, Phys. Fluids 14, 1476 (2002)], practical evaluation requires a 
Sonine polynomial approximation, and here it is mainly illustrated in the case 
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of the shear viscosity. Second, to check the accuracy of the Chapman-Enskog 
results, the Boltzmann equation is numerically solved by means of the Direct 
Simulation Monte Carlo (DSMC) method. The simulation is performed for 
a system under uniform shear flow, using the Gaussian thermostat to control 
inelastic cooling. The comparison shows an excellent agreement between the- 
ory and simulation over a wide range of values of the restitution coefficients 

and the parameters of the mixture (masses, concentrations, and sizes). 
PACS number(s): 45.70.Mg, 05.20.Dd, 51.10.+y 
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I. INTRODUCTION 



The macroscopic behavior of rapid granular flows can be described through hydrody- 
namic equations accounting for dissipation among the interacting particles. A basis for the 
derivation of the hydrodynamic equations and explicit expressions for the transport coeffi- 
cients appearing in them is provided by the corresponding Boltzmann kinetic theory in the 
low-density regime. In the simplest model the grains are taken to be smooth hard spheres 
with inelastic collisions. Assuming the existence of a normal solution for sufficiently long 
space and time scales, the Chapman-Enskog method [1] can be applied to get the velocity 
distribution function in terms of the hydrodynamic fields and their spatial gradients. This 
method must be conveniently adapted to inelastic collisions due to the new time depen- 
dence of temperature resulting from collisional energy loss. In the case of a monocomponent 
gas, the Navier-Stokes transport coefficients have been obtained in terms of the restitu- 
tion coefficient without limitation on the degree of inelasticity [2-4]. This analysis for a 
monocomponent system has been also extended to dense gases in the context of the Enskog 
equation [5]. 

Similar studies for multicomponent granular gases are more scarce and most of them 
limited to the asymptotically weak dissipation limit [6]. Although these studies permit in 
principle different temperatures for each species, they usually assume energy equipartition 
and so the partial temperatures Tj are made equal to the mixture temperature T. However, 
some recent results obtained in molecular dynamics simulations [7] as well as in real exper- 
iments of vibrated mixtures in three [8] and two [9] dimensions clearly show the breakdown 
of energy equipartition. A more recent kinetic theory calculation which takes into account 
temperature differences has been carried out by Garzo and Dufty [10]. They solved the set 
of Boltzmann coupled equations of the binary mixture by means of the Chapman-Enskog 
expansion for states near the local homogeneous cooling state. The mass, heat, and momen- 
tum fluxes were determined to first order in the gradients of the hydrodynamic fields and the 
associated transport coefficients were explicitly identified. As in the case of elastic collisions 
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[1], these transport coefficients verify a set of coupled linear integral equations which are 
solved approximately by using the leading terms in a Sonine polynomial expansion. The 
results derived by Garzo and Dufty [10] provide a description of hydrodynamics in binary 
granular mixtures valid a priori over the broadest parameter range and not limited to the 
quasielastic regime. In particular, the consequences of the temperature differences on the 
transport coefficients were shown to be quite significant. 

In the case of molecular fluid mixtures, it is known that the leading order truncation 
is quite accurate, except for extreme mass ratios (e.g., electron-proton systems). Much 
less is known in the case of inelastic collisions, although some comparisons with computer 
simulations for homogeneous states indicate that the accuracy is similar to that for elastic 
collisions [11,12]. The objective here is to compare the kinetic theory predictions for the 
shear viscosity with those obtained from a numerical solution of the Boltzmann equation 
by means of the Direct Simulation Monte Carlo (DSMC) method [13]. Specifically, the 
simulations are performed for a granular mixture undergoing uniform shear flow (USF), 
namely, a macroscopic state characterized by constant partial densities rii (i = 1, 2), uniform 
temperature T and a linear flow velocity profile u ijX = u x = ay, a being the constant shear 
rate. In a molecular fluid under USF, the temperature increases in time due to viscous 
heating. As a consequence, the average collision frequency v (which is proportional to 
T 1 / 2 for hard spheres) increases with time and the reduced shear rate a* = a/is (which 
is the relevant nonequilibrium parameter of the problem) tends to zero in the long time 
limit. This implies that for sufficiently long times the system reaches a regime described 
by linear hydrodynamics and the Navier-Stokes shear viscosity can be identified [14,15]. 
For granular fluids, the inelasticity of collisions introduces an energy sink in the balance 
equation for the temperature. Thus, the relationship between the temperature and the 
shear viscosity is not as simple as for molecular fluids since there is a competition between 
viscous heating and collisional cooling. However, if the granular fluid is externally excited 
by an external energy source that exactly compensates for the collisional energy loss, the 
viscous heating effect is still able to heat the system (as in the elastic case) and one can 
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identify the shear viscosity in the limit of small shear rate relative to the collision frequency 
v (i.e., a* — > 0). Although there are several choices for the external driving force, here we 
consider an external thermostat proportional to the peculiar velocity. This thermostat has 
been frequently used in nonequilibrium molecular dynamics simulations of molecular fluids 
[16] and has the advantage that, in the absence of shear, it does not affect the dynamics of 
the system at all since it is formally equivalent to a rescaling of the velocities [7] . 

The motivation of our study is twofold. First, in light of some doubts about the validity 
of a hydrodynamic description for granular flow, the comparison with simulation allows us 
to test the Chapman-Enskog solution obtained by assuming the existence of a normal or 
hydrodynamic regime. Since the parameter space here is quite large the tests of the theory 
and concepts are quite stringent. Second, as said above, we can also assess the degree 
of reliability of the approximate solution (first Sonine polynomial approximation) to the 
resulting integral equation over a wide range of the parameter space. With respect to the 
driving external force used in our analysis, we do not claim that it is the most suited one to 
model any real experiment. However, it has the advantage that it can be incorporated into 
kinetic theory and computer simulations very easily and it allows to check the assumptions 
of the Chapman-Enskog method. 

The plan of the paper is as follows. In Section II, we review the Boltzmann equation 
and associated macroscopic conservation laws in the presence of the Gaussian thermostat. 
The Chapman-Enskog method is applied in Section III to get all the transport coefficients 
of the mixture, with special emphasis in the shear viscosity coefficient. The details of the 
derivation are displayed in Appendix A. Section IV deals with the application of the DSMC 
method of the Boltzmann equation to USF with thermostat. The Chapman-Enskog and 
simulation results are compared in Section V at the level of the shear viscosity showing a 
good agreement. We close the paper in Section VI with a discussion of the results presented. 
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II. THE BOLTZMANN EQUATION AND TRANSPORT COEFFICIENTS 



We consider a binary mixture of smooth hard spheres of masses m,\ and m 2 , and diameters 
<7i and o 2 . The inelasticity of collisions among all pairs is characterized by three independent 
constant coefficients of normal restitution an, a 22 , and a± 2 = a 2 \, where is the restitution 
coefficient for collisions between particles of species i and j. Due to the intrinsic dissipative 
character of collisions, an energy supply is requested to fluidize a granular gas. For simplicity, 
here the fluidization is driven by the action of a non-conservative external force, frequently 
referred to as the Gaussian thermostat. In this case, the mixture is heated by an "anti-drag" 
force, linear in the peculiar velocity V and chosen to exactly compensate for collisional 
cooling. As said in the Introduction, this deterministic thermostat has been widely used in 
computer simulations of molecular fluids [16]. Under these conditions and in the low-density 
regime, the distribution functions /j(r, v; t) {i = 1, 2) for the two species are determined from 
the set of nonlinear Boltzmann equations 



where the constant £ is taken to be the same for each species [7,12,17]. Here, Vi = v x — u, 
u being the flow velocity. The Boltzmann collision operator Jij[vi\fi, fj] describing the 
scattering of pairs of particles is 



where = (<7j + <jj) /2, <r is a unit vector along their line of centers, 6 is the Heaviside 
step function, and g i2 = Vi — v 2 . The primes on the velocities denote the initial values 
{v^, v 2 } that lead to {v 1; v 2 } following a binary collision: 



(d t + Vl • v) h + i$ A • (Vx/i) = £ Jy [vxl/M, /,-(*)] , 



(1) 





x a^fiir, vi, t)fj(r, v 2 , t) - fi(r, Vi, t)fj(r, v 2 , t) 



(2) 





where /i^ = rrii/ (m ; + rrij). 



6 



The relevant hydrodynamic fields are the number densities n«, the flow velocity u, and 
the "granular" temperature T. They are defined in terms of moments of the distributions 

fi as 



rii 



j dvfi(v) , pu = J2j dvmivfi(v) , (4) 



nT=p = Y,fdv^V 2 f t (v), (5) 

where n = n\ + n 2 is the total number density, p = m x ni + m 2 n 2 is the total mass density, 
and p is the hydrostatic pressure. At a kinetic level, it is convenient to introduce the kinetic 
temperatures Tj for each species defined as 

\n lTl = j d^V'U (6) 

The collision operators conserve the particle number of each species and the total mo- 
mentum but the total energy is not conserved: 

J dvJij[v\fi, fj] = 0, Yl j dvrmvJijivlfi, f 3 ] = 0, (7) 

J2 I dvimiV 2 Jy[v|/i, /,] = -fnTC, (8) 

where ( is identified as the cooling rate due to inelastic collisions among all species. The 
macroscopic balance equations follow from the Boltzmann equation (1) and Eqs. (7) and 
(8). They are given by 

Dt-rii + n;V • u + — — = , (9) 

mi 

Au + p" 1 VP = 0, (10) 



n ; rrii on 

In the above equations, D t = d t + u ■ V is the material derivative, 



ji = m t J rfvV/,(v) (12) 
is the mass flux for species % relative to the local flow, 

P = E / dvmtVVfiiv) (13) 

i 

is the total pressure tensor, and 

q = W dv\ mi V 2 \ fi(v) (14) 

i J 

is the total heat flux. 

The energy balance equation (11) shows that the existence of a driving with the choice 
£ = C compensates for the cooling effect due to the inelasticity of collisions. In that case, 
the macroscopic balance equations look like those of a conventional mixture with elastic 
collisions. Nevertheless, the transport coefficients entering in the constitutive equations 
are in general different from those of a gas of elastic particles. Furthermore, for systems 
with elastic collisions, the specific set of gradients contributing to each flux is restricted by 
fluid symmetry, time reversal invariance (Onsager relations), and the form of the entropy 
porduction [18]. In the case of inelastic collisions only fluid symmetry applies and so there 
is more flexibility in representing the fluxes and identifying the corresponding transport 
coefficients. It follows from fluid symmetry that the pressure tensor has the same form to 
first order in the gradients as for the monocomponent gas. In the case of heat and mass 
fluxes, several different (but equivalent) choices of hydrodynamic fields can be used and some 
care is required in comparing transport coefficients in the different representations. Here, 
as done in the unforced case [10], we take the gradients of the mole fraction x\ = ni/n, the 
pressure p, the temperature T and the flow velocity u as the relevant ones. Thus, in this 
representation, the phenomenological constitutive relations for the fluxes in the low-density 
regime have the forms [19] 

m^ DVxi _ E to'VT, Ja = - Jl , (15) 

P p T 

q = -T 2 D"Vx 1 - LVp - AVT, (16) 



Paf3 = P5 af 3 - V yJflUa + V a M/3 - -5 af3 V • uj . (17) 

The transport coefficients are the diffusion coefficient D, the thermal diffusion coefficient D', 
the pressure diffusion coefficient D p , the Dufour coefficient D", the thermal conductivity A, 
the pressure energy coefficient L, and the shear viscosity 77. The Chapman-Enskog method 
[1] generalized to inelastic collisions allows one to get explicit expressions for these transport 
coefficients as functions of the restitution coefficients and the parameters of the mixture. 



III. SHEAR VISCOSITY OF A HEATED GRANULAR MIXTURE 

The Chapman-Enskog method assumes the existence of a normal solution in which all 
space and time dependence of the distribution function occurs through a functional depen- 
dence on the hydrodynamic fields 

fi(r, vi,i) = /itvilx^r,*),^,*),^,*),^^*)]. (18) 

This functional dependence can be made local in space and time by means of an expansion 
in gradients of the fields. Thus, we write fi as a series expansion in a formal parameter e 
measuring the nonuniformity of the system, 

/ i = /i 0) +6/i 1) +6 2 /f ) + -" , (19) 

where each factor of e means an implicit gradient of a hydrodynamic field. The local reference 
state fj ^ is chosen to give the same first moments as the exact distribution The time 
derivatives of the fields are also expanded as d t = + ed^ + ■ ■ ■. The coefficients of the 
time derivative expansion are identified from the balance equations (9)— (11) after expanding 
the fluxes, and the cooling rate ( in a similar series as (19). More details on the Chapman- 
Enskog method adapted to inelastic collisions can be found in Ref. [10]. Now, the main 
difference with respect to the free cooling case [10] is that the sink term in the energy 
equation is zero (when one takes £ = (), so that the terms coming from the time derivative 
vanish. 
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In the zeroth order, obeys the kinetic equation 



^•(Vif 0) ) = E^ 0) .ft (20) 

3 

where use has been made of the choice = Here, is determined by Eq. (8) to 
zeroth order. With this choice, Eq. (20) is identical to the the one obtained in the unforced 
case [20], and there is an exact correspondence between the homogeneous cooling state and 
this type of driven steady state. This is one of the advantages of the Gaussian thermostat. 
Dimensional analysis requires that //°^(V) must be of the form 

ri 0) (V)=n t v^(V/v ) (21) 

where 



„„ = 2T '^± (22) 
V mim 2 

is a thermal velocity defined in terms of the temperature T of the mixture. So far, the exact 
form of $j has not been found, although a good approximation for thermal velocities can be 
obtained from an expansion in Sonine polynomials [20] . In the leading order, $ j is given by 



9i 



3/2 



-OiV* 



IT 



1 + I (V* 4 - 56 t V* 2 + 



(23) 



where V* = V/v , &i = (/x,^) -1 , and 7, = Tj/T. The coefficients q (which measure 
the deviation of $j from the reference Maxwellian) are determined consistenly from the 
Boltzmann equation. The approximation (23) provides detailed predictions for the cooling 
rate the temperature ratio T 1 /T 2 and the cumulants q as functions of the mass ratio, size 
ratio, composition and restitution coefficients. Recently, the accuracy of this approximate 
solution has been confirmed by Monte Carlo simulation of the Boltzmann equation over a 
wide range of the parameter space [11] . 

The analysis to first order in e is similar to the one worked out in Ref. [10] for the free 
cooling case. Some details of the derivation of the transport coefficients as well as their 
final expressions are given in Appendix A. Here, given that the theoretical predictions for 
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the shear viscosity coefficient rj will be compared with those obtained from Monte Carlo 
simulations, we focus our attention on the explicit final expression for rj. According to Eq. 
(A25), the shear viscosity can be written as 



nT 
V = — V , 



(24) 



where 



y/nnal 2 v , 



(25) 



is an effective collision frequency and the reduced shear viscosity coefficient r/* is 



2 J* 



(26) 



with 



» _ 72(722 ~ C) - TiTjj - _ 7i(tii-C 

"1,1 ~ A ' "2,1 — A 



72721 



(27) 



Here, C* = C (0) M 



A = 7l72 [C 2 - C(rn + r* 2 ) + T n T 22 - t* 12 t* 21 



(28) 



and the dimensionless quantities t*- are given by [10] 



16 
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5V2 W/ 



1-4(1 -an) 



64 



+-^2i(l + [60f Vi20 2 - A*2i0i)(0i + d 2 y 1/2 

to L 

+l^ 2 (e l + e 2 f 2 (3 - a 12 ) + + o 2 r 1/2 



c 2 26 2 (l2fi 21 + 9/i 12 - 10) -9,(5- 6/^1) - f/x 2 i(3 - a 12 )(#i + 2 ) 



16 



# 2 ) 5 / 2 



(29) 



A*21 



'12 



. X2 _^(l + 6d 2 2 (^ 12 9 2 - ^iX^i + ^ 2 )" 1/2 

10 /ii2 L 



+^2i0 2 2 {9 1 + ^ 2 ) 1/2 (3 - a 12 ) - S^" 1 ^! + 9 2 ) 



-1/2 



ci 2^(10 - 12/i 12 - 9/i 2 i) + ^(5 - 6/i 12 ) - §//2i(3 - ai 2 )(0i + #2) 



16 



(^i + ^ 2 ) 5 / 2 



(30) 



11 



The corresponding expressions for r 2 * 2 and t 21 can be inferred from Eqs. (29) and (30) by 
interchanging 1 <-> 2. 

Equations (24)-(30) provide the explicit expression for the shear viscosity rj of a heated 
granular binary mixture in the first Sonine approximation. It is apparent that the reduced 
viscosity rf presents a complex nonlinear dependence on the restitution coefficients an, 
a 22 , and au and the parameters of the mixture mi/m 2 , <J\jo 2 , and n\jn 2 . The quality 
of the expression for 77* will be later assessed by comparison with Monte Carlo simulations 
in the uniform shear flow (USF) problem with a thermostat. Before studying the general 
dependence of rf on the parameter space, it is instructive to consider some special limit 
cases. In the elastic limit, an = a 22 = ^12 = 1, C* — 0> li — 1> and c\ = c 2 = 0, so that the 
expression (26) becomes 

„ _ XjRi + x^R 2 + X1X2-R12 
^ rrf 5i + + x 1 x 2 S 12 ' 



where 



= - + it> 2 = - + — , (32) 

3 5' 35^ 



8 v^2~ 
15 + T" 



°1 V -1/2 -1 , / ff 2 \ 2 -1/2 _l 

A% Vl2 + A*12 ^21 

(5\ 2 ) 



(33) 



s ' ^ 5^ l£J A ' 52 = wf * fei • (34) 

S -l + i(^) 2 (S) 2 ^" 1/2 ' <35 » 
Here, /x = /i 12/^21 = rn\jm 2 is the mass ratio. Equation (31) agrees with the results 
obtained in the first Sonine approximation to the coefficient of viscosity of a molecular gas- 
mixture of hard spheres [21]. In the case of mechanically equivalent particles (mi = m 2 , 
an = "22 = "12 = ol, 0i = <7 2 = cr), 7i = 1, C* = (2/3)(l - « 2 )(1 + 3c/32), and 

_ 32(l-q)(l-2a 2 ) 
Cl - C2 - C -81-17« + 30^(l-a) (36) 
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In this case, one gets 



rj 



15 

T 



((l + a) 



(2 + a) + — (33a -39) 



(37) 



This expression coincides with the one recently obtained for a heated granular monocompo- 
nent gas [4]. All this shows the self-consistency of the present description. 



IV. MONTE CARLO SIMULATION FOR UNIFORM SHEAR FLOW WITH 

THERMOSTAT 

The USF is a nonequilibrium state characterized by constant partial densities n i; a linear 
velocity profile u = Uj = a • r, where the elements of the tensor a are a« = a5k x Se y , a being 
the constant shear rate. In addition, the granular temperature T and the pressure tensor P 
are uniform, while the mass and heat fluxes vanish by symmetry reasons. This special state 
is generated by Lees-Edwards boundary conditions [22] which are simple boundary periodic 
boundary conditions in the local Lagrangian frame R = r — a • rt and V = v — a ■ r. In 
terms of these variables the velocity distribution functions are uniform [23] 

/i(r,v;0 = /i(V;*). (38) 

In the case of elastic collisions (( = 0) and in the absence of a thermostatting force, the 
energy balance equation (11) yields the heating equation 

d t T = ~^P xy . (39) 

Since the granular temperature T increases in time, so does the collision frequency u{t) 
according to Eq. (25). As a consequence, the reduced shear rate a*(t) = a/u(t) (which 
is the relevant uniformity parameter) monotonically decreases with time and the system 
asymptotically tends towards that of (local) equilibrium. This implies that for sufficiently 
long times (which means here a* < 1) the system reaches a regime described by linear 
hydrodynamics and the Navier-Stokes shear viscosity rj can be identified as [14,15] 
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where P* y = P xy /nT. This route has been shown to be quite efficient to measure the 
Navier-Stokes shear viscosity coefficient for dilute [24] and dense [15] gases. 

For a granular mixture, unless a thermostat is introduced, the energy balance equation 
(11) leads to a steady state when the viscous heating effect is exactly balanced by the colli- 
sional cooling [25]. However, if the granular mixture is excited by the Gaussian thermostat 



that exactly compensates for the collisional energy loss, the viscous heating still heats the 
system and Eq. (39) remains valid. Consequently, the linear relationship (40) allows one 
to determine the shear viscosity coefficient in the long time limit. Recently, this idea has 
been used to measure the shear viscosity of a heated granular monocomponent gas [4] . The 
comparison with kinetic theory shows an excellent agreement over a wide range of values of 
the restitution coefficient. It must be noted that here 77 represents the shear viscosity of an 
excited granular mixture and thus it does not necessarily coincide with the Navier-Stokes 
shear viscosity obtained in the unforced case [10]. As a matter of fact, the results obtained 
in Sec. Ill indicate that the transport properties are affected by the Gaussian thermostat 
and the expression (26) for the (reduced) shear viscosity differs from the one derived in the 
freely cooling case [10]. The use of thermostats to control collisional cooling in undriven 
systems is quite common [12,26]. Usually, the motivation is to produce a steady state while 
here is to remove the steady state in favor of one whose dynamics determines the viscosity. 

The Boltzmann equation for a mixture of inelastic hard spheres under USF and subject 
to the external Gaussian force (41) reads: 



The second term on the left-hand side represents an inertial force of the form F" 1 = — m^a ■ V, 
while the third term represents the thermostat force F* h given by Eq. (41). Thus, in this 




(41) 



dtU - *Vy W h + 2^ • W.) = E Ja [WMJAt)} 



(42) 
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frame, the system is in a homogenous state subjected to the action of the (total) force 
F- n + F- h . We have numerically solved Eq. (42) by means of the Direct Simulation Monte 
Carlo (DSMC) method [13]. This method was devised to mimic the processes involved in the 
Boltzmann collision term and its extension to deal with inelastic collisions is straightforward. 
In addition, since the USF is spatially homogeneous in the Lagrangian frame, the simulation 
method is easy to carry out and only the (peculiar) velocities of the particles need to be 
stored. The restriction to this homogeneous state prevents us from studying the possible 
formation of particle clusters (microstructure). 

Technical details of the DSMC method and its application to the USF state can be 
found in Refs. [11] and [25]. In our simulations we have typically taken a total number of 
particles N = 10 5 , a number of replicas Af = 5, and a time step At = 3 x I0~ 3 l n /v 01 , where 
£u = (V2nnicr'fj is the mean free path for collisions 1 — 1 and t> i — \j^T/mi. 

At given values of the shear rate a, the restitution coefficients and the parameters of 
the mixture, the system is initially prepared in a local equilibrium state with a temperature 
T(0) = T such that the initial value of the reduced shear rate is = a/u{T ). As the 
system evolves, we monitor the time evolution of the reduced shear rate a*(t) = a/v{T(t)) 
and the reduced xy element of the pressure tensor P* y (t) = P xy (t)/nT(t). We observe that 
in all the cases, after a transient period, the ratio rf = —P* y /a* reaches a constant value 
that is independent of the shear rate and time. This allows us to measure the corresponding 
shear viscosity coefficient r\ as 

,W = =gV, (43) 

where the dimensionless shear viscosity rf is independent of time but depends on dissipation 
and the parameters of the mixture (masses, sizes and concentrations). 

The theoretical prediction for rf can be obtained from the Chapman-Enskog solution to 
Eq. (42) up to first order in the shear rate a. In Appendix B it is easily proved that the 
first order solution to (42) leads to the same expression for the shear viscosity as the one 
obtained in Sec. Ill from the general Chapman-Enskog method specialized to USF. Thus, in 
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the first Sonine approximation, the theoretical prediction of r]* is given by Eqs. (26)-(30). 

Before analyzing the dependence of the dimensionless shear viscosity coefficient rf on the 
parameters of the problem, it is instructive to test the consistency of the simulation method 
in the limit a* — > (which corresponds here to tu >> 1) . For long times and for given values 
of a, mi/m 2 , 01/02 , and ni/n 2 , the reduced viscosity r/* must reach a value independent 
of the inital preparation of the system. In Fig. 1, we plot the shear-rate dependent viscosi- 
ties rj*(a*) measured in the simulation, relative to its Navier-Stokes value i]^ given by the 
Boltzmann theory [Eqs. (26)-(30)] for three different choices of the initial shear rate a^: 0.2, 
0.3, and 0.4. Here, the restitution coefficient is a = 0.9, the mass ratio is m\jmi = 4, the 
concentration ratio is ni/n 2 = 1 and the size ratio is 01/02 = 3. After a transient regime 
of a few mean free times, we observe that the curves corresponding to the three different 
initial conditions practically coincide. This means that a hydrodynamic regime independent 
of the inital conditions has been achieved. In addition, for very small values of a* 2 , the ratio 
rf{a*)/rfe fluctuate around 1 showing that in this regime the viscosity coefficient measured 
in the simulation is consistent with the value obtained from the Boltzmann kinetic theory. 
The same behavior has been found for other values of the restitution coefficient as well as of 
the parameters of the mixture. Notice that the limit a* — > is strictly unattainable in the 
USF because it requires an infinite amount of time. Also, the signal-to-noise ratio decreases 
in that limit so that the fluctuations increase. 



V. COMPARISON BETWEEN THEORY AND SIMULATION 

Once the consistency of the simulation method has been tested, we focus our attention on 
the study of transport properties in the Navier-Stokes regime. In this Section we compare 
the predictions of the Sonine approximation with the results obtained from the DSMC 
method. A complete presentation of the results is complex due to the high dimensionality 
of the parameter space: {an, CK22, CK12, mi/m,2, 01/02, ^1/^2}- For the sake of concreteness, 
henceforth we will assume that the spheres are made of the same material, i.e., an = a 22 = 
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«i2 = a. This reduces the parameter space to four quantities. 

Apart from the shear viscosity coefficient, another interesting quantity at this level of 
description is the temperature ratio Xi/T 2 . This ratio measures the breakdown of the energy 
equipartition. The analysis of the temperature differences has been a subject of growing 
interest in the past few years among both theorists [7,11,12,20] and experimentalists [8,9]. 
As was previously found from the Boltzmann kinetic theory [20], except for mechanically 
equivalent particles, the partial temperatures Tj are different. For the sake of illustration, 
Fig. 2 shows the dependence of the temperature ratio on the size ratio oxjoi for an equimolar 
mixture (ni/n 2 = 1) and three different values of the restitution coefficient a = 0.9, 0.8, and 
0.7. We consider a binary mixture of constant density and so, mi/m 2 = (o-i/a 2 ) 3 . We 
observe that for large size ratios the temperature differences are quite important, even for 
moderate dissipation. It is also apparent that an excellent agreement between the theory 
(given by the first Sonine correction) and Monte Carlo simulations (symbols) is found over 
the entire range of values of size and mass ratios considered. 

Next, we explore the influence of dissipation on the reduced shear viscosity r)*(a) for 
different values of the mass ratio, the size ratio, and the concentration ratio. Three different 
values of the (common) restitution coefficient are considered: a = 0.9,0.8, and 0.7. In 
Fig. 3, we plot the ratio rf{a)/rf{l) versus the mass ratio mi/m 2 for o\joi = n\jn-i = 1. 
Here, 77* (1) refers to the elastic value for the shear viscosity coefficient. Again, the symbols 
represent the simulation data while the lines refer to the theoretical results obtained from the 
Boltzmann equation in the first Sonine approximation. We see that in general the deviation 
of 7)* (a) from its functional form for elastic collisions is quite important. This tendency 
becomes more significant as the mass disparity increases. The agreement between the first 
Sonine approximation and simulation is seen to be in general excellent. This agreement 
is similar to the one previously found in the monocomponent case [4]. At a quantitative 
level, the discrepancies between theory and simulation tend to increase as the restitution 
coefficient decreases, although these differences are quite small (say, for instance, around 2% 
at a — 0.7 in the disparate mass case mi/m 2 = 10). 
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The influence of the size ratio on the shear viscosity is shown in Fig. 4 for mi/m 2 = 4 
and rii/n 2 = 1. We observe again a strong dependence of the shear viscosity on dissipation. 
However, for a given value of a, the influence of <7i/<72 on r]* is weaker than the one found 
before in Fig. 3 for the mass ratio. The agreement for both a = 0.9 and a = 0.8 is quite good, 
except for the largest size ratio at a — 0.8. These discrepancies become more significant as 
the dissipation increases (say a = 0.7), especially for mixtures of particles of very different 
sizes. Finally, Fig. 5 shows the dependence of rj* (a) / rj* (1) on the concentration ratio for 
mi/ 1712 = 4 and 0\joi = 1. We observe that both the theory and simulation predicts a 
very weak influence of composition on the shear viscosity. With respect to the influence of 
dissipation, the trends are similar to those of Figs. 3 and 4: the main effect of inelasticity in 
collisions is to enhance the momentum transport with respect to the case of elastic collisions. 
The agreement now between theory and simulation is very good, even for disparate values of 
the concentration ratio and/or strong dissipation. Therefore, according to the comparison 
carried out in Figs. 3, 4, and 5, we can conclude that the agreement extends over a wide range 
values of the restitution coefficient, indicating the reliability of the first Sonine approximation 
for describing granular flows beyond the quasielastic limit. 



VI. DISCUSSION 

Although the utility of a hydrodynamic description for granular media under rapid flow 
conditions has been recognized for many years, its domain of validity as well as the forms of 
the transport coefficients remain a topic of interest and controversy. In this context, there 
are some doubts about the possibility of going from a kinetic theory to a hydrodynamic level 
of description by using a Chapman-Enskog expansion around the homogenous cooling state. 
Given that the search for exact solutions of the Boltzmann equation is far beyond the present 
perspectives, an alternative to get some insight into the above question is to numerically solve 
the kinetic equation and compare these results with the corresponding solution obtained by 
assuming the validity of a hydrodynamic description. In this paper, we have performed such 
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a comparison at the level of the shear viscosity coefficient of a heated granular mixture. The 
system is heated by the action of a thermostatting external force which exactly compensates 
for cooling effects associated with the inelasticity of collisions. Although some previous works 
[4,27,28] have compared kinetic theory predictions for transport coefficients with computer 
simulations in the case of a monocomponent gas, studies for multicomponent granular gases 
are more scarce. Very recently, a seemingly similar analysis for the shear viscosity r\ of a 
dense mixture has been given in Ref . [29] . Nevertheless, the above kinetic theory only holds 
for nearly elastic particles and the expression of rj in the first Sonine approximation coincides 
with the one obtained in the elastic case. 

As a first step in our issue, in Sec. Ill we have derived the general hydrodynamic equations 
of a heated binary mixture of smooth inelastic spheres from the Boltzmann kinetic equation 
by using the Chapman-Enskog method. The corresponding transport coefficients have been 
expressed in terms of the solution to integral equations, which are then solved approximately 
(first Sonine polynomial approximation) just as in the case of elastic collisions. The explicit 
expressions for the set of relevant transport coefficients {D, D p , D', D", L, A, 77} are displayed 
in Appendix A. In contrast to previous works [6,29], our results are not limited a priori to 
weak inelasticity and they take into account the effect of the temperature differences on 
the transport coefficients. On the other hand, the results obtained here for the transport 
coefficients slightly differ from those obtained in the freely cooling case [10], showing that 
in general the introduction of a thermostat affects the transport properties of the gas [23]. 
The Chapman-Enskog results obtained for the mixture have been then specialized to the 
hydrodynamic state of transverse shear. In this state, the shear viscosity coefficient r\ is the 
relevant transport coefficient of the problem. The explicit form of rj is given by Eqs. (26)- 
(30) in terms of the restitution coefficients and the parameters of the mixture (masses, 
diameters, and concentrations). 

To test the assumptions of the Chapman-Enskog method and the approximate Sonine 
solution to the resulting integral equation, the DSMC method has been used to solve the 
Boltzmann equation in the uniform shear flow state. In the absence of a thermostat, in 
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a granular fluid there is a competition between two opposite effects: viscous heating and 
collisional cooling. In that case, when both effects exactly cancel each other, a steady state is 
reached after a transient period. In this steady state, due to the coupling between dissipation 
and the shear rate, the system is far away from the Navier-Stokes regime, except when a — > 1 
[25]. However, if the external thermostat is adjusted to compensate for the energy lost in 
collisions, the shearing work still heats the system. As a consequence, as the system evolves, 
the reduced shear rate a*(t) goes to zero and the system achieves a regime described by 
linear hydrodynamics. In this regime, the Navier-Stokes shear viscosity coefficient can be 
measured from simulations. In this paper, the thermostat is used to remove the steady state 
in favor of a time-dependent state whose dynamics allows one to get the Navier-Stokes shear 
viscosity just as for the case of elastic collisions [15,24]. 

The dependence of the viscosity rj on the full parameter space has been explored. Specifi- 
cally, the parameter space over which our solution has been verified is the mass ratio mi/m 2 , 
the concentration ratio ni/n 2 , the ratio of diameters ai/a 2 , and the (common) restitution 
coefficient a = an = an = a 22 = a i2 - The theory and simulation clearly show how in 
general, the influence of dissipation on momentum transport is quite important since there 
is a relevant dependence of the viscosity 77(a) on the restitution coefficient a. At a given 
value of the restitution coefficient, the dependence of r)(a)/r)(l) on the mass ratio is more 
significant than the one found on the composition and diameters. This feature has been also 
found for the temperature ratio in the experiments recently carried out in vibrated mixtures 
[8,9], although experimental confirmation of the trends observed here for the viscosity is still 
lacking. With respect to the accuracy of the theory predictions we see that, in general, the 
Chapman- Enskog results in the first Sonine approximation exhibit an excellent agreement 
with the simulation data. This supports the idea that the Sonine polynomial approximation 
for granular fluids has an accuracy comparable to that for elastic collisions. Exceptions to 
this agreement are extreme mass ratios and strong dissipation. These discrepancies are due 
basically to the approximations introduced in applying the Chapman-Enskog method, and 
more specifically in using the first Sonine approximation. 
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One of the main limitations of the results obtained here from the Boltzmann equation is 
its restriction to the low-density regime. In this regime, the collisional transfer contributions 
to the fluxes are negligible and only their kinetic contributions are taken into account. Pos- 
sible extension in both aspects, theory and simulation, of the present simple hydrodynamic 
state to higher densities can be carried out in the context of the revised Enskog theory. In 
this case, many of the phenomena appearing in dense granular fluids (such as spontaneous 
formation of dense clusters surrounded by regions of low-density [30]) could be studied. On 
the other hand, although the comparison performed here has been made undergoing uniform 
shear flow without paying attention to the possible formation of particle clusters [31], our 
Chapman- Enskog results apply for general inhomogeneous situations. The only restriction is 
that they provide the irreversible parts of the mass, heat, and momentum fluxes to leading 
order in the spatial gradients of the hydrodynamic fields. In this context, the results derived 
in this paper can be used to analyze the behavior of granular mixtures in a lot of physical 
situations. Thus, for instance, the knowledge of the complete hydrodynamic equations for 
a binary mixture allows one to say whether the mixture hydrodynamics is more or less un- 
stable to long-wavelength perturbations than that of the one-component case, and what are 
the mechanisms involved in phenomena very often observed in nature and experiments such 
as phase separation or segregation. We hope that the present results give some insight into 
the understanding of these interesting and complex problems. 
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APPENDIX A: CHAPMAN-ENSKOG EXPANSION 

In this Appendix, the expressions of the transport coefficients for a heated granular 
mixture are obtained. The derivation follows similar steps as those made in Ref. [10] in the 

21 



free cooling case. Here, we will use the same notation as in Ref. [10]. In the first order, the 
distribution function ff verifies the kinetic equation 



(d(" + v.v)/M, 



where D*jp = df + u • V, and 



(U f (o)j 



(Al) 



(A2) 



Miff' = -^[./i ,/j ]■ 



(A3) 



In these equations, it is understood that i ^ j and use has been made of the fact that 
d\ 0) T = and the results jf } = q(°) = C (1) = 0. The last equality follows from the fact 
that the cooling rate is a scalar, and so should be proportional to V • u. However, as 
shown later, there is no contribution to ff proportional to the divergence of the flow field 
so that = by symmetry. This property is special to the low density Boltzmann kinetic 
theory and such terms occur at higher densities [5]. The macroscopic balance equations to 
first order are 

D { px 1 = 0, -D^ Inp = ^A (1) InT = -V • u, T> t (1) u = -p _1 Vp. (A4) 
5 2 

Use of these in (Al) yields 

(A + ^C (0) ^ • ff ] + Miff = A t • Vxx + Bj • Vp + Cj • VT + A^V*^, (A5) 
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Q(V) = - 
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The solutions to Eqs. (A5) are of the form 



fW = Ai ■ Van +Bi-Vp + Ci-VT + V iiaP V a up 



(A10) 



The coefficients A%,Bi,Ci, and T> ij0l p are functions of the peculiar velocity V and the hy- 
drodynamic fields. The cooling rate depends on space through its dependence on x±, p, 
and T. The integral equations for the unknowns are easily identified as coefficients of the 
independent gradients in (A10). The result is 



Ai 



'c- + -c<®—-v 
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Bi 

Note that, in contrast to what happens in the free cooling case [10], here each one of the 
quantities Ai, Bi, C iy and T> ia p obey closed integral equations. The solution to Eq. (All) 
provides the expression for the transport coefficients. In the case of the mass flux ji, these 
coefficients are identified as 



D = 



P 



3m 2 n 



J dvV-Ai, 



mip 



D p = -_^|rfvV -B ± , 
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3p 



The transport coefficients for the heat flux are 



(A12) 



(A13) 
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\hj dv\m z V 2 V- Bi, 
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A = 4E / dv\ mi V 2 V -C t . (A17) 

i=l 1 

Finally, the shear viscosity is given by 

1 2 r 

r? = --E / dvrmVaVpVirf. (A18) 
iu i=1 J 

Accurate aproximations to the solutions to the integral equations for (Ai, Bi, Ci, 2\ Q| g) 
may be obtained using low order truncation of expansions in a series of Sonine polynomials. 
In the case of the mass flux, we consider the leading Sonine approximation (lowest degree 
polynomial) 

{Ai,Bi,Ci} -> AmVKiA.i,^}, / ilA f(V) = n l K/2vrT J ) 3 / 2 exp(-m i \/ 2 /2T i ), (A19) 

where a i: i = —(m 1 m2n/pn i T i )D, b iy i = —(p/priiT^Dp, and c i:1 = -(p/TrnT^D 1 . The 
coefficients a^i, b^i, and c^i are determined by multiplying the three first equations of Eq. 
(All) by m^V and integrating over the velocity. The result is 

a M =-(, D 4c<°f (^1-,^, (A20) 

ci,i = 0. (A22) 

Here, the collision frequency z/ D is given by Eq. (73) of Ref. [10]. 

In the case of the pressure tensor, the leading Sonine approximation for the function 

A,a/3 — * fi,Mdi t lRi t a/3, Ri,a/3 = mi(V a Vj3 — \V 2 5 a p). (A23) 

The shear viscosity coefficient is given by 

V = -nT 2 (xiTiXi + ^27 2 2 4,i) • (A24) 
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The coefficients d iy \ can be determined by multiplying the fourth equation of Eq. (All) by 
Ri^afi and integrating over the velocity to get the coupled set of equations 

\ 



/ r n - C (0) r 12 ^ ' 



K T21 r 22 - c (0) 



di,i 

y ^2,1 j 
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The collision frequencies 77,- = t*z/, where t*- are given by Eqs. (29) and (30). From Eq. 
(A25) one easily gets the expression (24) for the shear viscosity given in the main text. 

The calculations for the heat flux are similar to those previously made for the other 
fluxes. As in the unforced case, this requires going to the second Sonine approximation. In 
this case, the transport coefficients defining the heat flux (16) are given by 
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The coefficients 2 , fej 2 and Cj 2 obey the equations 
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The corresponding expressions of the elements X 2 , Y" 2 and Z 2 can be deduced from Eqs. 
(A32), (A33) and (A34), respectively, by interchanging 1^2 and setting D — > D, D p — > 
—Dp and £>' — > — £>'. The frequencies Uij and the collision integrals appearing in Eqs. (A31), 
(A32), (A33) and (A34) were explicitly evaluated in the Appendix D of Ref. [10]. Thus, the 
transport coefficients D", L, and A are completely determined. 



APPENDIX B: FIRST ORDER SOLUTION TO THE USF 

In this Appendix we get the solution to Eq. (42) in the first order in the shear rate a. 
The normal solution to Eq. (42) is provided by the Chapman-Enskog method, i.e., a solution 
given as a power series in a: 

U = ff ) + f ( ^ + --- (Bl) 

The zeroth-order solution ff^ verifies Eq. (20) and it corresponds to the homogeneous cooling 
state distribution in the local Lagrangian frame. Its first Sonine approximation is given by 
Eq. (23). Inserting the expansion (Bl) into Eq. (20) leads to the following integral equation 
for if: 

MP - aV y Jrif + ±C< 0) ^ • (V// 1 )) + i C (1) 7^ • (V/D = - M,rr. (B2) 
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where the operators and Mi are defined by Eqs. (A2) and (A3), respectively. Since f- ' 
depends on time only through the temperature, Eq. (39) implies that dtf} 1 ^ = 0(a 2 ) and so 
the first term on the left hand side of Eq. (Bl) vanishes in the first order. Further, = 
by symmetry because V • u = in the USF. Taking into acount the above properties, Eq. 
(B2) reduces to 

(C + ±C<»'^ ■ v) /« + Mjf = aV y A/<«). (B3) 

This integral equation is identical to Eq. (A5) when one particularizes the latter one to 
USF. Therefore, the expression for the shear viscosity obtained from (B3) is given by Eqs. 
(24)-(30). 
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FIGURES 

FIG. 1. Plot of the ratio rf{a*)/rf^ as function of a* for a = 0.9 in the case m\/m2 = 4, 
ni/n-i = 1, and oi/o<i = 3 for three different values of the initial shear rate a^: Uq = 0.2,0.3, 
and 0.4. Here, i]^ refers to the Navier- Stokes shear viscosity value given by the first Sonine 
approximation to the Boltzmann equation. 

FIG. 2. Plot of the temperature ratio T1/T2 as a function of the size ratio a\/(T2 = (mi/m^) 1 / 3 
for n\jn2 = 1 and three different values of the restitution coefficient a: (a) a = 0.9 (circles), (b) 
a = 0.8 (squares), and (c) a = 0.7 (triangles). The lines are the theoretical predictions and the 
symbols refer to the results obtained from the DSMC method. 

FIG. 3. Plot of the ratio rj*(a)/r)*(l) as a function of the mass ratio m\/m,2 for 
<7i/<72 = n i/n2 = 1 and three different values of the restitution coefficient a: (a) a = 0.9 (circles), 
(b) a = 0.8 (squares), and (c) a = 0.7 (triangles). The lines are the theoretical predictions and 
the symbols refer to the results obtained from the DSMC method. 

FIG. 4. Plot of the ratio rj* (a) / rf (1) as a function of the size ratio <J\jo2 for m\/m2 = 4, 
n\/n2 = 1 and three different values of the restitution coefficient a: (a) a = 0.9 (circles), (b) 
a = 0.8 (squares), and (c) a = 0.7 (triangles). The lines are the theoretical predictions and the 
symbols refer to the results obtained from the DSMC method. 

FIG. 5. Plot of the ratio rf (a) / rf* (1) as a function of the concentration ratio n\jri2 for 
mi/ m2 = 4, <ti/<72 = 1 and three different values of the restitution coefficient a: (a) a = 0.9 
(circles), (b) a = 0.8 (squares), and (c) a = 0.7 (triangles). The lines are the theoretical predic- 
tions and the symbols refer to the results obtained from the DSMC method. 
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